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Abstract 

The loop gas approach to lattice field theory provides an alternative, geometrical 
description in terms of fluctuating loops. Statistical ensembles of random loops 
can be efficiently generated by Monte Carlo simulations using the worm update 
algorithm. In this paper, concepts from percolation theory and the theory of self- 
avoiding random walks are used to describe estimators of physical observables 
that utilize the nature of the worm algorithm. The fractal structure of the random 
loops as well as their scaling properties are studied. To support this approach, 
the 0(1) loop model, or high-temperature series expansion of the Ising model, is 
simulated on a honeycomb lattice, with its known exact results providing valuable 
benchmarks. 

Keywords: loop gas, Monte Carlo, worm update algorithm, fractal structure, 
critical properties, duality 



1. Introduction 

Representing the hopping of particles from one lattice site to the next, the 
strong-coupling expansion in relativistic quantum field theories formulated on a 
spacetime lattice provides an alternative approach to numerically simulating lat- 
tice field theories in terms of world lines. The standard approach, which is rooted 
in the functional integral approach to field quantization, involves estimating ob- 
servables (expressed in terms of the fields) by sampling a representative set of field 
configurations. New configurations are typically generated by means of a Monte 
Carlo technique which uses importance sampling, with each field configuration 
weighted according to the probability that it occurs. In contrast, the approach 
based on the strong-coupling, or hopping expansion, which is closely connected 
to Feynman's spacetime approach to quantum theory yj], involves linelike objects. 
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Figure 1 : Bond-shifting algorithm for generating a new world line configuration on a cubic lattice. 
Lattice sites visited by the walk are marked by full circles and the updated plaquettes are shaded. 

Physical observables are in this geometrical approach no longer estimated by sam- 
pling an ensemble of field configurations, but by sampling a grand canonical en- 
semble of (mostly closed) world lines, known as a loop gas, instead. The weight 
of a given world line configuration is typically determined by the total length of 
the paths, the number of intersections, and the number of loops contained in the 
tangle. 

In statistical physics, the strong-coupling expansion is known as the high- 
temperature series expansion yO. Lattice field theories studied in this context 
are typically spin models, such as the O(iV) spin model, whose representation in 
terms of high-temperature (HT) graphs is known as a loop model. 

A first numerical study of loop gases formulated on the lattice was carried 
out by Berg and Foerster yfl. New world line configurations were generated by 
a bond-shifting Monte Carlo update algorithm as follows. A randomly chosen 
bond of the existing configuration is shifted perpendicular to itself by one lattice 
spacing in any of the 2(d— 1) directions of the hypercubic lattice. During the shift, 
each of the endpoints of the moving link erases or draws a bond in the chosen 
perpendicular direction, depending on whether the link is occupied or not, as in 
Fig. ED The new configuration is accepted or rejected according to the Metropolis 
algorithm. 

At about the same time, Dasgupta and Halperin fl, following a suggestion 
by Helfrich and Miiller ^\ that the HT graphs of the O(N) lattice model simul- 
taneously describe a loop gas of sterically interacting physical lines, simulated a 
gas of directed loops on a cubic lattice. New loop configurations were generated 
in this study by inserting an elementary loop, or plaquette, of random orientation 
according to the Metropolis algorithm. 

Although these and related early loop gas update algorithms |@, 0, work 
fine in the disordered phase away from the critical point, they all, being based on 
local updates, suffer from pronounced critical slowing down. That is, consecutive 



configurations are highly correlated close to the critical point, and simulations on 
larger lattices become increasingly unfeasible in this region. 

About a decade ago, Prokof'ev and Svistunov [9] have introduced a Monte 
Carlo update algorithm that, although based on local updates, does away with 
critical slowing down almost completely. The so-called worm algorithm generates 
loop configurations, not by inserting plaquettes, but through the motion of the end 
points of an open world line — the "head" and "tail" of a "worm". An additional 
loop is generated in this scheme when the head bites the tail, or through a "back 
bite" where the head erases a piece (bond) of its own body and thereby leaves 
behind a detached loop and a shortened open chain. 

Besides this outstanding technical advantage, the worm algorithm has the ad- 
ditional advantage in the context of statistical physics that the complete set of 
standard critical exponents can be determined at a stroke. This set is known to 
split into two, viz- the thermal and the magnetic exponents. While the thermal 
exponents, such as the specific heat exponent a, pertain to closed paths, the mag- 
netic exponents, such as the magnetic susceptibility exponent 7, pertain to open 
paths in the geometrical approach. Using a plaquette update, one is restricted to 
the topology of the initial configuration. If that starting configuration consists of 
just closed paths, a plaquette update algorithm will subsequently also generate 
only loop configurations. Open paths, needed to determine the magnetic expo- 
nents, must be sampled in such a scheme by putting in an open path connecting 
two fixed endpoints from the start. A plaquette update will then change the loops 
fluctuating in the background and will also change the form of the open path, but 
it will leave the endpoints of the path untouched. Since, in principle, all possible 
end-to-end distances are needed to determine the magnetic exponents, a plaque- 
tte update is impracticable to achieve this. By the nature of the worm algorithm, 
which features an open path between loop updates, these data are generated on the 
fly in this scheme. More specifically, the open paths directly sample the spin-spin, 
or two-point, correlation function. 

In this paper, which extends previous work by two of us on the subject II 1 Ol . 1 1 IN , 
we describe estimators of physical observables that naturally arise in a loop gas 
and that allow determining the standard critical exponents. Our approach, put 
forward in Sec. |2l amalgamates concepts from percolation theory — the paradigm 
of a geometrical phase transition — and the theory of self-avoiding random walks. 
We relate this geometrical approach to phase transitions in terms of fluctuating 
paths to the more familiar field theory approach by considering the O(N) sym- 
metric 4 theory in Sec. [3] To support our arguments, the second part of the paper 
is devoted to Monte Carlo simulations of the two-dimensional 0(1) loop model 



3 



using the worm update algorithm. This model serves as a prototype with its var- 
ious exact results providing a yardstick for our Monte Carlo results and also for 
the feasibility of our approach. Section |4] specifies the model we have simulated, 
introduces the specific implementation of the worm update algorithm used, and 
gives details of the simulations. Our results are presented in Sec. [5J We finish 
with a discussion and outlook. 



2. Loop gases 

We are concerned with lattice field theories close to the critical point where 
they undergo a continuous phase transition. Their equivalent loop gas representa- 
tion can be conveniently characterized by the average number £ n of closed paths, 
or polygons, of n steps per unit volume. Close to the critical point K r , the so- 
called loop, or loop length distribution takes asymptotically a form [10] similar 
to the cluster distribution near the percolation threshold known from percolation 



theory UJ], 



~ n-^-V 6 ™, 9 oc (K - K c ) 1/a . (1) 



Here, 9 is the line tension (in suitable units), K is the tuning parameter, and d 
denotes the dimension of space (in the case of classical theories) or spacetime (in 
the case of quantum theories). When the line tension is finite, the Boltzmann fac- 
tor in the distribution (OQ) exponentially suppresses long loops. Upon approaching 
the critical point, 9 vanishes at a rate determined by the exponent a. At K c , loops 
proliferate for they can now grow without energy penalty. The remaining factor in 
the loop distribution is an entropy factor, giving a measure of the number of ways 
a polygon of n steps can be embedded in the lattice. It is characterized by the 
fractal dimension D of the paths at the critical point. The entropy factor decreases 
with increasing n. 

A standard definition of the fractal dimension is through the asymptotic be- 
havior of the average square radius of gyration (R^) of chains of n steps as 

(Rl) ~ (2) 



where 



k,k'=l k=l 



with x ik the position vector of the chain after k steps and 

1 n 



x 

n 

k=l 
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the center of mass of the chain (which can be closed or open) of n steps. Here 
and in the following, lattice sites are labeled by the index i. The radius of gyration 
gives a measure of the distance covered by the path. Another standard definition 
is through the average square end-to-end distance (Rf) of open chains of n steps, 

(R 2 c ) = ((x ln -x l0 ) 2 )~ n 2 / D , (5) 

where Xi denotes the starting point of the chain. For the two-dimensional O(N) 
model, which for — 2 < N < 2 undergoes a continuous phase transition, the 



fractal dimension D of the HT graphs II 1711 corresponds to the renormalization 



group eigenvalue y 2 of the two-leg operator in the spin representation of the model 



The natural length scale in quantum field theory is the correlation length £. 
The critical exponent u, characterizing the divergence of this length scale when 
the criticalpoint is approached, £ ~ \K — K c \~ u , is related to the fractal dimension 



through 



u = 1/aD. (6) 



This expression, which assumes the same form as in percolation theory 111211 . gen- 
eralizes a celebrated result due to de Gennes [14] for self-avoiding random walks 
(SAWs), which corresponds to the limit N — > of the O(iV) spin model. In that 
case, cr = l, but in general a takes different values, see Tabled] below. 

As is known from the theory of SAWs, closed paths alone yield only the ther- 
mal exponents of the universality class defined by the 0(iV — > 0) model. To obtain 
also the magnetic exponents, and thereby the complete set of standard exponents, 
the total number 

Z n = ^ ^ Xj) (7) 

j 

of SAWs of n steps starting at Xi and ending at an arbitrary site Xj is needed 
in addition. Because of translational symmetry, z n does not depend on Xj, and 
z n (xi,Xj) only depends (up to lattice artifacts) on the end-to-end distance r = 
\xi — Xj\, i.e., z n (xi,Xj) = z n {r). The ratio of z n (xi,Xj) and z n defines the 
probability P n (xi, Xj) of finding a chain connecting X{ and Xj in n steps. As for 



SAWs Ill5h . we expect this distribution to scale for a general loop gas as 

P n {xi,x 3 ) = z n (x hXj )/z n ~ n^^Vir/n 1 ^), (8) 

with V a scaling function. That is, we assume that P n (xi,Xj) depends only on 
the ratio r/(R^} 1 ^ 2 . As an aside, the average square end-to-end distance © is the 
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second moment of this distribution. In continuum notation: 



(Rl) = n d / drr d - l r 2 P n (r) 
Jo 



(9) 



where Q d denotes the surface of a unit hypersphere embedded in d space dimen- 
sions, and P(r), being a probability, is normalized to unity 

POO 

l = Q d drr^Pnir). (10) 
Jo 

In addition to the scaling ([8]), we also assume the number z n to scale as 

z n K n ~ vf' D e- 6n (11) 

with a universal exponent i3 that characterizes, as do the rest of the critical expo- 
nents, the universality class. For the 0(N) model, it depends, in addition to the 
dimensionality d, solely on N. Since the number of possible rooted open chains 
with no constraint on their endpoint increases with the number n of steps, d is 
expected to be positive. This is in contrast to closed chains, where the corre- 
sponding factor in Eq. (Q~|) decreases with increasing n, reflecting that it becomes 
increasingly more difficult for chains to close the longer they are. 

The fractal dimension D together with the exponents a and i? determine the 
standard critical exponents of the theory. As for SAWs, the relevant scaling rela- 
tions can be derived by writing the correlation function G{x il Xj) as a sum over 
all possible chains of arbitrary many steps joining the endpoints: 

G(Xi, Xj) = ^ Z n(^i, Xj)K n . (12) 

I! 

As before, G(xi, Xj) = G(r) because of translational invariance. When evaluated 
at the critical point, where the correlation function depends algebraically on the 
end-to-end distance, G(xi, Xj) ~ l/r d ~ 2+v , this gives 

rj = 2- D-d. (13) 



Given the exact values for 77 II16I1 and the fractal dimension D of the HT graphs 
[fTvh . $ can be determined exactly for the two-dimensional O(N) model, see Ta- 
bled] Through the exact enumeration and analysis of the number z n of SAWs on 
a square lattice up to length 7 1 , the expected value {) / D = ^for N = Q has been 



established to high precision Ill8h 
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Table 1: Critical exponents of the two-dimensional critical O(iV) spin models, with N = 
—2, — 1, 0, 1, 2, oo, respectively, together with the fractal dimension D of the HT graphs as well as 
the two exponents a and 

7 



Model 


N 


Gaussian 


-2 




-1 


SAW 





Ising 


1 


XY 


2 


Spherical 


oo 



37 
32 
43 
32 

7 

1 

OO 
OO 



V 


V 


D 


a 


d 





i 


5 


8 


3 


2 


4 


5 


4 


3 


5 


13 


16 


11 


20 
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10 


13 


20 


5 
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11 


21 
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21 
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11 


8 


3 


1 


8 


11 


8 


1 
i 


OO 


3 
2 





1 
1 





OO 


2 









The relation © can, incidentally, be derived by using the second-moment 
definition of the correlation length £, 

C 2 _ J™ dr r*- 1 r 2 G(r) 

S TOO , J 1 \ l^^V 



/ °°drr^ 1 C7( 



in continuum notation. Finally, using the definition of the susceptibility \-> X — 

^2 j G(xi, Xj), which diverges as x ~ 1-^" — -^c|~ 7 , we find 

7 = (D + $)/aD. (15) 



This relation generalizes one originally due to des Cloizeaux 111911 for SAWs for 
which a = 1. The explicit expressions for u, r), and 7 satisfy Fisher's scaling 
relation, j/u = 2 — rj. Note that only the combinations D + $ and aD enter the 
scaling relations between the various critical exponents. 

We next consider the limit Xj — > Xi of z n (xi, Xj). Following standard practice 
in the theory of SAWs 1I20I1 . we define this limit of vanishing end-to-end distance 
as the number of chains z n (xi,Xi ± afi) = z n (a) of n steps returning to a site 
Xi ± afi adjacent to the starting point Xi. Here, ±fi is a unit vector in any (pos- 
itive as well as negative) direction on the lattice (see Fig. [2]), and a is the lattice 
spacing. That is, z n (a) rather than z n (0) is taken when closing open chains, with 
the lattice spacing a serving as an ultraviolet cutoff. In continuum quantum field 
theory, the limit x' — > x corresponds to putting two fields at the same point. Such 
composite operators usually need special care and require a separate renormaliza- 
tion independent of that of the constituting operators. The number of chains z n (a) 
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Figure 2: A SAW on a square lattice returning to a site adjacent to its starting point x. 



is related to the loop distribution (OQ) through 

1 



n 



z n (a)K n . (16) 



Since a polygon can be traced out starting at any lattice site along the chain, the 
factor 1/n is included to avoid double counting. Note that the loop distribution is 
a density being defined per lattice site and that z n (a) refers to rooted closed chains 



all starting at the same lattice site As first shown by McKenzie and Moore 112111 
for SAWs, consistency of Eq. (fT6l) with z n (a) = z n P n (a) and Eq. (OQ) requires that 
the scaling function V(t) must vanish for t — > and behave for small argument t 
as 

V{t) ~ f (17) 

with an exponent determined by the asymptotic behavior (fTTT) of the number z n 
of open chains at the critical point. With this identification, Eq. (fT3l) becomes the 
relation first proposed by Prokof 'ev and Svistunov in Ref. ll2~2ll . Together with the 



relation © proposed in Ref. [HQ, Eq. (PT3T) allows expressing the standard critical 
exponents in terms of the fractal structure of open and closed paths and the rate 
I /a at which the line tension vanishes upon approaching the critical point. As 
already mentioned in the Introduction, a major advantage of the worm update al- 
gorithm is that it features both open and closed paths because this makes possible 
to determine all these exponents at a stroke. 

3. 1 4 1 Lattice field theory 

To make connection with field theory, we consider as an example the O(N) 
symmetric 4 theory formulated on a hypercubic lattice in d Euclidean spacetime 
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dimensions. The theory is specified by the (Euclidean) lattice action 



S = ad E { i E + - + ^V(^) + f ^) I (is) 

with lattice spacing a. The real scalar field </?( x i)> which is defined on the lattice 
sites Xi of the spacetime box, has N components </? = y? a = y? 2 , . . . , <^ Ar ). 
As before, the index i labels the lattice sites, the sum J2i stands for a sum over 
all lattice sites, and the index a — 1, 2, . . . , N labels the field components. More- 
over, ip 4 = (<£> • y?) 2 , where the dot product implies a summation over the field 
components: ip ■ ip = 5^=1 < £' C V Q! - Lattice coordinates, representing discretized 
spacetime, are specified by )i, with /t denoting the unit 

vector pointing in the (positive) ^-direction. Moreover, m 2 is the bare mass pa- 
rameter squared, and g is the bare coupling constant of the self-interaction term. 
In the world line picture, this four-leg operator corresponds to intersections where 
two lines cross. The renormalization group eigenvalue y 4 of the four-leg operator 
corresponds to the fractal dimension D x of these intersections. Numerically, this 
fractal dimension can be determined through finite-size scaling by measuring the 
average number, or "mass" M x , of these intersections which scales at the critical 
point as 

M x (L)~L Dx (19) 

with the linear size L of the lattice. For the critical two-dimensional O(N) model, 
the four-leg operator is irrelevant for — 2 < N < 2, i.e., y 4 = D x < 0, and 



it becomes marginal for N = 2 fl 1 3f1 . Because intersections are irrelevant (or 
marginal) there, loops at the O(iV) critical point are frequently referred to as dilute 
loops. 

In the continuum limit, where the lattice spacing tends to zero, a — > 0, the 
lattice action (fT8l reduces to the standard form 

S = jd d x^ [dM*)} 2 + ^<P 2 (x) + f ¥> 4 (*)} , (20) 

where (p(x) stands for the field defined in continuous spacetime. 

The partition function Z of the lattice theory obtains by carrying out the sum, 
or integral over the spin variable at each site of the lattice: 

Z = Tre~ s , (21) 
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with 



Tr = J] J d N <p( Xi ). (22) 



This amounts to summing, or integrating over all possible spin configurations, 
each weighted by the Boltzmann factor e~ 5 (in natural units). In the continuum 
limit a — > 0, this defines the functional measure J Dip, and the partition function 
becomes 

Z= ( Dpe~ s . (23) 



For numerical simulations, a more convenient form of the lattice action is ob- 
tained by casting Eq. (IT8T ) in terms of dimensionless fields and parameters defined 



through MM 



a d "V(^) = 2tf# (24) 

a*~ d g = 6 A (25) 

m V = 1 ~l XN - 2d, (26) 
K 

with K > 0. The action then takes the form of an O(iV) spin model 

s = ■ <k + + A E (ti - N ) 2 ■ w 

(*,*') « « 

The sum Yltn') extends over all nearest neighbor pairs. In terms of these new 
dimensionless variables, the action is independent of the lattice spacing a. The 
partition function Z can now be written as 

Z= I D/^expfi^^A (28) 

V <v'} / 

with the on-site measure 

J DM) = J JJd^ e^- x ^- N T. (29) 



In the limit A — > oo, the lattice field theory reduces to the standard O(iV) spin 



model, with a "spin" variable of fixed length, (ff = N, located at each site of 
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the spacetime lattice. The remaining factor in the on-site measure (|29l) becomes 
trivial in this limit and can be ignored. The normalization is chosen such that 



1 N & = 1, Jd N ^ = N. (30) 

Instead of considering the conventional Boltzmann weight factor, often a sim- 
plified representative of the O(N) universality class is studied, obtained by trun- 



cating that factor I124ll : 



Z= / J] J] (1 + (31) 

* (i,i'> 

The second product is restricted to nearest neighbor pairs. The main difference 
with the original spin model is that in the truncated model, links cannot be multi- 
ply occupied. The weight carried by a configuration is positive for \K\ < l/|iV|. 
By universality, the truncated model is expected to still belong to the O(N) uni- 
versality class. Note that for N — 1, where 0j = ±1, the full Boltzmann factor 
can be exactly written in the truncated form by the identity 

ffiMi' = cosh(/3) [1 + tanh(/3)0;0;,] oc 1 + Kfcfo (32) 

with K = tanh(/3). The prefactor cosh(/3) is immaterial and can be ignored as 
far as critical phenomena are concerned. The worm algorithm [0] was originally 
designed to simulate the HT representation of the theory (l28l) with the full Boltz- 
mann factor included so that links can be multiply occupied. However, as already 
suggested by its inventors [@], the algorithm can be readily adapted to simulate 
the truncated model ( f3Tb without multiple occupied links. 

The scaling part of the logarithm of Z reads expressed in terms of the loop 
distribution 

lnZ/y~^4, (33) 

n 

with V the volume. The result © immediately follows from the hyperscaling 
argument that In Z/V ~ 

In (continuum) quantum field theory, the two-point correlation function G(x, x') 
is given in the symmetric phase by the average of a product of two ip fields at x 
and x', respectively: 

G{x,x') = {cp(x)-^{x')}. (34) 
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Its algebraic behavior at the critical point is in this context parameterized as 

G(x,x') ~ \/ r 2d f with 

d v = \{d-2 + -q) = \{d-D- ■&) (35) 

denoting the anomalous scaling dimension of the cp field. The limit x' — > x of the 
correlation function G(x, x') is conventionally defined through a "mass insertion" 



as da] 



d_ 

dm 

By Eqs. d33l and (|26|) it then follows that for a loop gas 



G(0) = (<f 2 {x)) oc — ^ In Z. (36) 



(<^(x)> ~ (K c - Kf/^Y.ntn, (37) 



or 



& 2 (x))~l/&\ (38) 

with 

d 2 = d- - (39) 
z/ 

the standard expression for the scaling dimension of the composite operator (p 2 (x). 
In deriving this, use is made of the relation ©. 

Note that naively taking the limit x' — > x in Eq. (PT21) yields, after using 
Eq. (TTBI) . the result (1371) without prefactor. This is the world line counterpart of 
the observation that composite operators usually require a multiplicative renor- 
malization by themselves that cannot be expressed in terms of the renormalization 
factors of the constituting operators. Also the power-law decay £F7J) of the scaling 
function V(t) for t — > is related to this. Assuming that the scaling function 
remains finite in the limit t — > 0, one obtains from Eq. (fT2l) with the relation © 
the incorrect result 

G(0) ~ Yl z nPn(a)K n rsj 1/^, (incorrect) (40) 

n 

involving the anomalous dimension of if instead of (/? 2 . Only for noninteracting 
theories, the renormalization of composite operators can be expressed in terms of 
the renormalization of the constituting operators, and V(0) is nonzero. 

Noting that the right side of Eq. (136*1 ) physically denotes the internal energy, we 
conclude from Eq. (|37l ) that in the world line approach this quantity is determined 
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by the average number of bonds in closed graph configurations, i.e., by the average 
total loop length in configurations without open chains. 

In closing this section, we remark that the two combinations D + -d and aD, 
on which the standard critical exponents depend, determine the anomalous scaling 
dimension of the (p and ip 2 fields through Eqs. (1331) and (|39l ) with l/u — aD, 
respectively. 



4. Model and details of simulation 

4.1. Loop model 

To specify the model we have simulated, we start with the representation (I3TI ) 
of the O(iV) model. Expanding the product appearing there, one readily verifies 
that only terms with an even number of spins at each lattice site contribute to the 
partition function. A factor K(f>f(f)j (no summation over a) in such a term can 
be conveniently visualized by drawing a bond along the link of the underlying 
lattice connecting the nearest neighbor sites labeled by i and j. With each field, 
or spin, component a = 1, 2, . . . , N is associated a color, so that the bonds come 
in N colors. Terms contributing to Z then correspond to closed graphs made up 
of such bonds and of vertices connecting an even number of bonds. The partition 
function is obtained by adding all these contributions, i.e., by summing over all 
possible disconnected closed graph configurations, each carrying a certain weight. 

The spin-spin, or two-point, correlation function G{x il Xj) of the truncated 
model 

G{ Xi , = (4>i -<f>j)=z f n & ■ fa n t 1 + R fa' ■ fa') (4i) 

can be treated in a similar fashion as the partition function with the proviso that for 
terms in the expansion of the product in the numerator to contribute, the two sites 
labeled by i and j must house, in contrast to all other lattice sites, an odd number 
of spins. Graphically, such terms typically correspond to a set of disconnected 
closed graphs with an additional open graph connecting the two endpoints and 

Xj. 

The O(N) loop model is obtained by resolving each closed graph into a unique- 
ly defined set of possibly intersecting loops. This is done by providing instructions 
how vertices connecting more than two bonds are to be resolved. In principle, 
such "walking instructions" can be formulated on an arbitrary lattice in arbitrary 



dimensions [25]. However, the simplest way to deal with this issue is to consider 
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a honeycomb lattice, which has coordination number z = 3, so that closed graphs 
simply cannot intersect. A configuration Q of disconnected closed graphs then 
automatically decomposes into loops, and the partition function of the resulting 



loop gas assumes the form Il24h 



Z loop = Y,K h N\ (42) 

Q 

where b denotes the number of bonds and I the number of loops in the graph. 
Each bond in a graph configuration carries a weight K, while each loop carries a 
degeneracy factor N, for they can have any of the N colors. These factors play 
the role of bond and loop fugacities in the loop model. The number of bonds in 
a graph configuration increases with increasing bond fugacity K and vice versa. 
The critical point of the 0(N) loop model on a honeycomb lattice is exactly known 



to be given by lll6ll 



K c = [2 + (2 - iV) 1 / 2 ] 1/2 . (43) 

4.2. Update algorithm 

In the main simulations, we restricted ourselves to the Ising model (N = 
1) on a honeycomb lattice. From a loop gas perspective, this model defines a 
statistical ensemble of polygons built from bond variables b t which are defined on 
the links of the lattice. Reflecting the fact that only one color is present (N = 
1), the bond variables only take the values b\ — 1, when the bond labeled by / 
is set, or bi = 0, when it is not. Apart from considering loop configurations, 
we also consider configurations that have in addition a single chain connecting 
two endpoints, X{ and Xj say. Such configurations, which correspond to two spin 
insertions in the spin representation, contribute to the numerator Z(xi, Xj) of the 
spin-spin correlation function 

G(x llXj ) = ^p± (44) 

and are naturally generated by the worm algorithm that locally updates the bond 
configurations. Although polygons on a honeycomb lattice cannot intersect, an 
open chain can "back bite" or touch a polygon. Such configurations, where a 
chain endpoint connects three bonds, are allowed and must be included in the 
update scheme. For a general loop model, such configurations pose a problem, 
for they can lead to a change in the number of loops during the next bond update. 
Then to keep track of the number of loops, the open chain must be traced out 
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anew, making the update algorithm nonlocal and slowing it down considerably. 
The Ising model is special in that the loop fugacity is unity, so that configurations 
with lattice sites housing three bonds do not pose a problem, at least not when 
just updating and not measuring them (see below). For the 0(1) loop model, 
the updates involve Metropolis flips of single bonds where the value bi of the 
bond variable is replaced with 1 — 6;. During the Monte Carlo simulation, chain 
endpoints move and, thus, accumulate information about open chain properties, 
such as their end-to-end distance. As there is a finite probability for an open chain 
to close and form a polygon, the algorithm also acquires information about the 
loop gas. We adapted the original worm algorithm [91] as follows, see Refs. [26, 



2711 for related adaptations. 

For configurations containing, in addition to polygons, a single chain with 
end-to-end distance larger than one lattice spacing, the updating scheme proceeds 
by 

1 . randomly choosing either endpoint of the chain, 

2. randomly choosing any of the links attached to the chosen endpoint, 

3. updating the corresponding bond variable bi with a single-hit Metropolis flip 
proposal bi — > b\ = 1 — bi with acceptance probability 

P accept = min(l,^ 1 - 2fei ) (45) 

as can be inferred from the weight K b in the partition function (l42l) . assum- 
ing that < K < 1. The exponent \ — 2bi = ±1 denotes the difference in 
the number of bonds contained in the proposed and the existing configura- 
tions. It follows that a proposal to create a bond is accepted with probability 

-Paccopt = K (< 1), whereas a proposal to delete one is always accepted. 

These updates are simple and straightforward as long as the chain remains open. 
Once, however, the chain has an end-to-end distance of just one lattice spacing, the 
existing configuration can be turned into a loop gas configuration by a single bond 
flip. Such an update then connects two different sectors of the model. Namely, the 
sector with an open path, which samples the numerator Z(xi, Xj) of the correlation 
function (l44l) . and the loop sector, which samples the partition function Z. In their 
original work [9J], Prokof'ev and Svistunov introduced conditional probabilities, 
parameterized by < p < 1, for Monte Carlo moves between the two sectors. 
We in this work put this parameter to unity and, thus, always attempt to close 
such a chain by using the update scheme above with the Metropolis acceptance 
probability (l45l) . If the update is accepted, and the open chain turns into a polygon, 
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we proceed by randomly choosing one link among all links of the lattice. The 
bond variable on that link is then subjected to a Metropolis trial move with the 
acceptance probability (1451) . 

To check the correctness of this worm algorithm, we simulated the critical 
0(1) loop model on a small 5x5 square lattice with periodic boundary conditions, 
i.e., on a torus and measured the spin-spin correlation function. Table |2] summa- 
rizes our Monte Carlo results and compares them to the exact results, obtained by 
direct enumeration. The table shows complete agreement within statistical error 
bars. 

4.3. Lattices 

In our main study, the loop model (1421) is regularized on a two-dimensional 
honeycomb lattice. As remarked before, the coordination number of the honey- 
comb lattice is three and allows a unique decomposition of closed graphs into an 
ensemble of polygons. We constructed the honeycomb lattice from its dual, i.e., 
hexagonal, or triangular, lattice. The latter, which, in contrast to the former, is a 
Bravais lattice, is spanned by two vectors of equal length, making an angle of 60°. 
We have chosen the lattice spacing of the dual lattice to be unity, a& = 1. The lat- 
tice spacing of the honeycomb lattice is then fixed to be = 1/ v^3 = 0.5773 . . .. 
Euclidean distances on the honeycomb lattice are measured in units of ciq. The 
number of lattice sites on the dual lattice, i.e., the volume, is taken to be V/\ = L 2 , 
where L denotes the number of lattice sites in any of the two independent direc- 
tions on the hexagonal lattice. This parameter L features as the linear lattice size 
variable in all our further considerations, including our finite-size scaling analy- 
ses. Under the dual construction, the volume of the honeycomb lattice picks up 
a factor of two, so that Vq = 2V/\, while the number B of links is unchanged, 
B A = B = 3V. 

We constructed a honeycomb lattice that is compact with periodicity of 2L in 
three directions. Let A M (xj) denote the shift operation that connects a site Xj to 
its nearest neighbor in the /x = 1,2,3 direction. Then there exist three product 
operations, each involving 2L such shifts, that form the identity and map any site 
Xi onto itself, see Fig. [3] for an example with L = 16. 

4.4. Observables 

In our Monte Carlo simulations, we analyzed the two sectors of the model, 
i.e., configurations with and without an open chain, separately. We implemented 
a search algorithm that uniquely decomposes (disconnected) closed graphs into 
polygons along the links of the underlying honeycomb lattice. Each polygon is 
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Table 2: Precision check of our worm update algorithm for the critical Ising model on a square 
5x5 lattice with periodic boundary conditions. We simulated the two-point correlation function 
G(xi,Xi) by putting one spin at the origin, labeled by 1, and the other spin at the site labeled by 
the index i, which goes through the lattice in a typewriter fashion. The second column gives results 
from exact summations, while the third column summarizes our Monte Carlo data obtained with 
the worm algorithm. The fourth column shows our Monte Carlo data in units of the exact results, 
and the last column gives the bit- variable t that is unity if theory and simulation differ by less than 
one a, and is zero otherwise. 



% 




G{x ll x i ) 


G / Gcxact t 


1 


1.000000 


1 





1.0 1 


2 


0.768360 





768353(30) 


0.999991(39) 1 


3 


0.708394 





708385(23) 


0.999987(33) 1 


4 


0.708394 





708342(38) 


0.999927(54) 


5 


0.768360 





768354(32) 


0.999993(41) 1 


6 


0.768360 





768350(40) 


0.999987(52) 1 


7 


0.722100 





722082(38) 


0.999976(52) 1 


8 


0.695433 





695370(33) 


0.999910(48) 


9 


0.695433 





695422(39) 


0.999986(56) 1 


10 


0.722100 





722175(31) 


1.000103(42) 


11 


0.708394 





708360(35) 


0.999952(49) 1 


12 


0.695433 





695412(37) 


0.999970(53) 1 


13 


0.683390 





683478(26) 


1.000129(38) 


14 


0.683390 





683430(46) 


1.000060(67) 1 


15 


0.695433 





695342(43) 


0.999870(62) 


16 


0.708394 





708390(40) 


0.999994(57) 1 


17 


0.695433 





695401(44) 


0.999955(63) 1 


18 


0.683390 





683318(42) 


0.999895(61) 


19 


0.683390 





683426(31) 


1.000053(46) 


20 


0.695433 





695408(51) 


0.999965(73) 1 


21 


0.768360 





768363(30) 


1.000005(39) 1 


22 


0.722100 





722037(36) 


0.999912(49) 


23 


0.695433 





695449(44) 


1.000024(64) 1 


24 


0.695433 





695465(43) 


1.000046(61) 1 


25 


0.722100 





722112(44) 


1.000017(60) 1 
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Figure 3: Compact honeycomb lattice with periodicity in three directions. The three operations 
mapping the site Xi marked by a full circle onto itself through shift operations are indicated by 
the three polygons, each winding the lattice once. Note the set bond in the upper right corner, 
which by the periodicity of the lattice belongs to the polygon winding the lattice in the northwest 
direction. 



assigned a length parameter n, denoting the number of sites visited. While tracing 
the loops, we also record for each site x ik visited by a polygon, the direction to 
the next site x ik+1 in that polygon. This makes possible to determine whether a 
polygon winds the (periodic) lattice in any of the three possible directions. Specif- 
ically, we determine for each polygon its winding number, which is a topological 
invariant, telling how often it winds the (periodic) lattice in a given direction. Note 
that because the endpoints of the worm erase or draw bonds, the worm algorithm 
can change the winding number of a configuration. A nonzero winding number is 
the signal for loop percolation. As in percolation theory, such "infinite" polygons 
are usually excluded from measurements of variables not connected to percolation 
observables to facilitate finite-size scaling analyses. 

For observables analyzed on lattices of fixed size, an even stringent upper 
bound on chain lengths is required. To this end, we monitored during the Monte 
Carlo runs the length n of winding loops. The minimum loop length found in the 
time series, each involving 10 7 sweeps of the lattice, represents a natural upper 
bound on (open and also closed) chain lengths to be included in such analyses. 
The minimum values n at the critical point for the honeycomb lattice of several 
sizes are given in Table |3j The length n (L) can be equally well interpreted as 
the length of the largest loop that can be realized within a lattice of linear size L. 
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Table 3: Length no of the shortest loop winding the honeycomb lattice of linear extent L recorded 



during 


_0 7 sweeps of the lattice at the critical point. 


L 


32 64 96 128 160 192 


224 


256 


288 


320 


352 


n 


78 196 346 510 690 854 


1138 


1320 


1602 


1820 


2022 



Stated differently, n (L) indicates the loop length up to which the scaling law (OQ) 
applies, see Fig[4j Scaling implies that this length increases with the linear lattice 
size as uq{L) ~ LP . The results in Table [3] satisfy this scaling with D = y, as 
anticipated, see Tabled! 

We measure the loop distribution t n by compiling a histogram of loop lengths 
during long Monte Carlo runs 

in = jr Yl 5 w (46) 

m 

where m enumerates the polygons measured with J2 m = N denoting the total 
number of polygons measured, and n m is the length of the mth polygon. Figure |4] 
shows the results of such measurements at the critical point on the largest lattice 
considered, i.e., L = 352. Loops at all scales are observed. The bump at the end 
of the distribution followed by a rapid falloff is typical for such distributions mea- 
sured on a finite lattice with periodic boundary conditions. Finally, we determine 
for each polygon its center of mass © as well as the square radius of gyration ©. 

We next turn to the analysis of configurations containing an open chain in ad- 
dition to polygons. Unlike configurations without one, those with an open chain 
cannot always be uniquely decomposed, even on a honeycomb lattice. Consider, 
for example, an open chain with an endpoint housing three bonds. It is then not 
clear whether it represents a single, self-intersecting chain, or a chain and a sep- 
arate polygon which touch each other at the chain endpoint. The decomposition 
is unique on a honeycomb lattice only if both chain endpoints contain just one 
bond. To minimize the arbitrariness associated with tracing out non-unique open 
chains, we omitted from the measurements chains with both endpoints housing 
three bonds. This only concerns a small fraction of all chain configurations mea- 
sured of about two percent at the critical point. When only one endpoint of the 
open chain contains three bonds, we interpret the configuration as representing 
an open chain and a separate polygon. This is similar in spirit to the closing of 
open SAWs, see the argument leading to Eq. (fT6l) and Fig. [2l As for polygons, we 
determined the square radius of gyration R 2 of open chains as a function of the 
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1e+010 




Figure 4: Log-log plot of the loop distribution l n as a function of the loop length n on the hon- 
eycomb lattice of size L = 352 at the critical point. The arrow indicates the minimal length 
no = 2022 from Table[3] The straight line proportional to n -2 /- 0-1 with D = ^ is put through 
the data points by hand to show the expected behavior ((TJ. 
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chain length n, as well as their square end-to-end distance R% also as a function 
of n. 

4.5. Details of simulation 

The simulations are carried out on honeycomb lattices of linear extent rang- 
ing from L = 32 to L = 352 in steps of AL = 32. Each simulation consists 
of 10 7 sweeps, where a single sweep is defined as Vq = 2Va = 2L 2 local bond 
updates. An additional 10% of the sweeps is used for thermalization. The ac- 
cumulated computer time used for the simulations amounts to a few weeks on a 
single workstation. 



5. Simulation results 

To check our code, we make use of the celebrated Kramers-Wannier duality 



for the two-dimensional Ising model 112811 . This duality asserts that the loop gas, 



or HT, representation of the 0(N = 1) model on a two-dimensional lattice at the 
same time represents the model in the standard spin representation on the dual 
lattice. To picture a spin configuration on the dual lattice, imagine drawing bonds 
between any pair of nearest neighbor spins on that lattice which are in the same 
spin state. When no bond exists between two nearest neighbor spins, they are in 
different spin states. A HT bond in a given loop configuration can be interpreted 
as indicating a broken bond between the two nearest neighbor spins living on the 
dual lattice, on either side of the HT bond on the original lattice. That is, a HT 
bond indicates that the two corresponding spins on the dual lattice are in different 
spin states. Given this transcription, a loop can then be pictured as forming the 
boundary of a cluster of nearest neighbor spins on the dual lattice which are all 
in the same spin state. This implies that loop configurations containing just a 
few bonds correspond to ordered spin configurations on the dual lattice. More 
generally, under the dual map, the high-temperature phase of the loop model, 
in which, for sufficiently high temperatures, only a few small loops are present, 
maps onto the low-temperature phase of the spin model, where, for sufficiently 
low temperatures, large spin clusters can be found. The two temperatures can be 
related by noting that a HT bond carries a factor K = tanh(/3), while a nearest 
neighbor pair on the dual lattice of unlike spins on each side of the HT bond carries 



a Boltzmann weight exp(— 2/3), so that Il28ll 



K = e~ 2 K (47) 
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Figure 5: Probability distribution function (PDF) of the magnetization (M) measured at K c , i.e., 
at the critical point of the infinite lattice, on lattices of size L = 32. One data set, marked by x , is 
obtained using the Swendsen-Wang cluster algorithm on the hexagonal lattice, the other data set, 
marked by +, is obtained using the worm algorithm on the honeycomb lattice, and then transcribed 
to the hexagonal lattice. Both sets are seen to nicely blend, as expected by duality. 



Note that the dual map as described here is special to two dimensions and cannot 
be generalized to higher dimensions. 

Not any loop configuration can be resolved in a spin configuration on the dual 
lattice. Loop configurations containing, for example, a single loop winding the 
lattice once do not, given the periodic boundary conditions, translate into a spin 
configuration on the dual lattice. It is straightforward to see that when the winding 
number of a given loop configuration is even in any of the three directions, such a 
transcription is possible up to a factor Z(2) which is chosen at random. 

To demonstrate duality and also the correctness of our Monte Carlo simula- 
tions, we measured the magnetization M of the Ising model at the critical point, 
using the two representations. For the standard spin representation of the Ising 
model on the hexagonal lattice, we use the Swendsen-Wang cluster update Il29ll . 
while for the loop model, or HT representation, on the honeycomb lattice, we use 
the worm algorithm. Only those loop configurations are considered that can be 
mapped onto a spin configuration on the hexagonal lattice. Figure [5] attests that, 
as expected, the two distinct data sets nicely merge. 
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Figure 6: Average absolute value of the magnetization on a hexagonal lattice as a function of the 
inverse temperature f3. The curve gives the exact result d49b due to Potts. 



As a further illustration, we display in Fig.[6]the average of the absolute value 
of the magnetization \M\ as obtained from the dual map as a function of the in- 
verse temperature (3 of the Ising model on the hexagonal lattice introduced in 
Eq. (|47l) . The simulation itself was carried out on the L = 96 honeycomb lattice 
for inverse dual temperatures /3 > /3 C , where by Eqs. (|47l ) and ( |43l 



f3 c = ±ln3 = 0.27465... . 



The plotted curve corresponds to the exact calculation by Potts BOh 



M 8 (/3) 



1 - 



16e -12/3 



(l + 3e- 4 /3) (1 -e- 4 ^)' 



(48) 



(49) 



for the functional form of the magnetization of the Ising model as a function of 
the inverse temperature f3 on an infinite hexagonal lattice. The agreement between 
the theoretical curve and the data is seen to be excellent. 

As a final check on the correctness of our Monte Carlo simulations, we mea- 
sured the Binder parameter 



U, 



1 

1 - - 



(M 4 ) 
3(M2) 



2 • 



(50) 



23 



0.6128 

0.6126 

0.6124 

0.6122 

0.6120 h 

0.6118 

0.6116 

0.6114 

0.6112 

0.6110 



hexagonal 
honeycomb 



Yl- si' 



50 100 150 200 250 300 350 400 

L 

Figure 7: Binder parameter Ul as a function of the linear lattice size L. One data set is obtained 
using the Swendsen-Wang cluster algorithm on the hexagonal lattice, the other set is obtained using 
the worm algorithm on the honeycomb lattice, and then transcribed to the hexagonal lattice. The 
straight line indicates the high-precision result obtained in Ref. ll3lll by transfer matrix methods. 



which involves the second and fourth moments of the magnetization M, see Fig.fTJ 
The definition of this parameter is such that for a Gaussian theory it vanishes. 
The critical value of the Binder parameter of the Ising model on a hexagonal 
lattice with periodic boundary conditions is known extremely well from transfer 
matrix calculations pill , viz. Ul = 0.61182773(1). Having no scaling dimension, 
Ul does not change with lattice size in leading order. One-parameter fits to the 
data lead to the estimates U L = 0.611822(73) with xV DOF = 1.45 for the loop 
gas on the honeycomb lattice (transcribed to the hexagonal lattice), and U L = 
0.611815(17) with x 2 /dof = 1.01 for the standard spin representation on the 
hexagonal lattice. Both estimates, based on different representations of the Ising 
model and obtained using different update algorithms, are in excellent agreement 
with the high-precision result. 

We next turn to estimators of physical observables which exploit the nature of 
the worm update algorithm and which can be naturally measured in this scheme. 
As first estimator we introduce the binary variable that records whether a loop 
configuration can be mapped onto a spin configuration on the dual lattice, or not. 
If a map exists, this observable is assigned the value zero, else it is assigned the 
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value unity. The top panel in Fig. [8] shows the average Il(/3) of this observable 
as a function of the inverse dual temperature /3 introduced in Eq. (1471) on lattices 
of linear extent L = 32, 64, 96. For large (3, where mostly only a few small loops 
are present, loop configurations can typically be mapped onto spin configurations 
on the dual lattice and Il(p) is small, tending to zero in the limit (3 — > oo. When 
(3 is lowered, larger loops appear and eventually loops can be found that wind 
around the lattice. As mentioned above, when, for example, a single loop does 
so once, a transcription is impossible, and the value unity is recorded in the time 
series of measurements. This explains the increase of Il((3) with decreasing (3. 
In the limit /3 — > 0, where loops are abundant, this observable is seen to tend 
to an asymptotic value (very close to) |. In words, the ratio of the number of 
loop configurations with an odd winding number to those with an even winding 
number tends to | in the zero-temperature limit, where it is recalled that, although 
impossible for loop configurations with an odd winding number, configurations 
with an even winding number can be mapped onto a spin configuration on the 
dual lattice. The observable Il{(3) has no scaling dimension and plays a role 
similar to the Binder parameter. Finite-size scaling implies that it depends not on 
the inverse temperature (3 and the lattice size L independently, but only on the 
combination (/3 — (3 c )L l / v with v the correlation length exponent. The bottom 
panel in Fig. [81 which displays the same data, but now as a function of this scaling 
variable with v — 1, shows that finite-size scaling is satisfied. 

To quantify this statement, we measured at the critical point on lattices 

of different size, see inset of Fig. |9] A one-parameter fit to the data leads to the 
estimate I L (/3 C ) = 0.50024(21) with x 2 /dof = 0.851. That is, half of the critical 
configurations have an odd winding number. Since Il((3) depends only on the 
scaling variable (f3 — fi^L 1 ^, differentiation of Il(/3) with respect to (3 allows 
estimating 1/V, see Fig. [9] A two-parameter fit to the data gives as estimate for 
the slope \ jv = 1.0001(15) with x 2 /dof = 1.01, in excellent agreement with the 
expected value v = 1 . 

As mentioned above, the great virtue of the worm algorithm is that it also 
generates open chains. Figure HO] shows the distribution z n , introduced in Eq. ©, 
of chains of n steps with arbitrary end-to-end distance r measured at the critical 
point on a lattice of linear size L = 352. ByEq. (fl~2"l) and the definition of the 
magnetic susceptibility x in terms of the correlation function given below Eq. (fl~4l) . 
the sum ^ n z n gives x at the critical point. According to finite-size scaling, x{K c ) 
scales with lattice size L as 

X (K c ) = Y,Zn~L^. (51) 
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Figure 8: Top panel: Average Il{P) as a function of the inverse dual temperature f3 an lattices 
of linear extent L = 32, 64, 96. Bottom panel: Same data reploted as a function of the scaling 
variable (fi - ft^L 1 ^ with v = 1. In the limit /3 ->• 0, tends to f (straight line). 
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Figure 9: Inset: Binary observable Il(/3 c ) as a function of the linear lattice size L. Main panel: 
Derivative I' L ($ C ) of with respect to (3 evaluated at the critical point f3 c as a function of L. 



The inset of Fig. [TO] shows a log-log plot of the ratio L 2 j\ as a function of the 
lattice size. A linear two-parameter fit to the data gives the estimate 2 — j/u = 
0.2498(26) with x 2 / D0F = 1.87, in agreement with the exact result r\ = 2 - 

7/" = I- 

An important characteristic of the chains generated by the worm algorithm, 
whether closed or open, is their Hausdorff, or fractal, dimension D. Figure [TT] 
shows the average square radius of gyration (R 2 ) of closed and open chains as 
a function of n, as well as the average square end-to-end distance (R%) of open 
chains, also as a function of n. In obtaining an accurate estimate of the fractal 
dimension from these and corresponding data measured on lattices of different 
size, we face two restrictions. The first is that the scaling © only holds asymp- 
totically, i.e., for sufficiently large n. To lift this restriction somewhat, we include 
the leading correction to scaling in Eq. © by writing 

(4e> = <™ 2/D (l + h A . (52) 



As for SAWs 113211 . we expect this leading correction term to be analytic and in- 
versely proportional to n. The length of a chain is, on the other side, bounded 
by the number of lattice sites available on a finite lattice. Figure [10] shows that, 
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Figure 10: Distribution z n of open chains of arbitrary end-to-end distance at K c on a lattice of 
linear size L = 352 as a function of chain length n. Inset: Log-log plot of the inverse of the 
integrated observable z n — x divided by the volume L 2 , i.e., of L 2 jx, measured on lattices 
of different size, as a function of L. 



as a result, even at the critical point, the number z n of open chains of arbitrary 
end-to-end distance falls exponentially for large n on a finite lattice. To obtain an 
estimate for the infinite lattice, chains that exceed a maximum length no where 
they start to notice the finite extent of the lattice must be ignored. As measure 
of this maximum, we take the length of the shortest polygon winding the lattice 
recorded in the time series, see Table [3j To be on the safe side, the maximum 
value of n included in the data considered for fitting is chosen to be about 0.8 n , 
while the minimum value it taken to be n min = 100. Table |4] summarizes the esti- 
mates for the fractal dimension D obtained through three-parameter fits using the 
form <l52l ). The last line in the table gives the weighted average of the estimates 
for L > 200, which are all consistent with a constant, i.e., L-independent value. 
The final results are in excellent agreement with the exact value -D C xact = y pre- 
dicted by Saleur and Duplantier Il33n . Drawing on numerical work by Cambier 
and Nauenberg B4I1 . Vanderzande and Stella 113 511 provided early indirect support 
for this prediction. Direct numerical support was first provided by Dotsenko et 
al. Q36D, and more recently in Ref. IllOn using other, percolationlike estimators and 
a plaquette update. 
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Table 4: Estimates of the fractal dimension D obtained through three-parameter fits to the average 
square radius of gyration (i?g) of open chains (left) and to the average square end-to-end distance 
(right) obtained on lattices of size L. 



L 


D 


D / D cxac t 


X 2 /dof 


D 


D / -D exact 


X 2 /dof 


160 


1.485(84) 


1.080(61) 


0.424 


1.43(14) 


1.04(10) 


0.682 


192 


1.4213(62) 


1.033(19) 


0.493 


1.475(46) 


1.073(33) 


0.588 


224 


1.3822(96) 


1.0052(70) 


0.301 


1.381(16) 


1.004(11) 


0.288 


256 


1.3672(67) 


0.9943(49) 


0.499 


1.364(12) 


0.9924(90) 


0.453 


288 


1.3762(44) 


1.0008(32) 


0.867 
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Figure 11: Average square end-to-end distance (scaled by a factor of five for readability) of 
open chains, as well as the average square radius of gyration (R^) of open and closed chains, all 
shown as a function of their length n and measured at the critical point on a lattice of linear size 
L = 352. 
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ln(i) 

Figure 12: Log-log plot of the average length (n w ) of loops winding the honeycomb lattice of 
linear extent L. 

Polygons that wind the lattice have been excluded from the above measure- 
ments. However, as in percolation theory, the fractal dimension of polygons (clus- 
ters) at the critical point can also be determined by exclusively focussing on such 
winding polygons (percolating clusters). Because these polygons are long, no cor- 
rections to scaling as in Eq. (l52l) need to be included. Figure [T2l shows the average 
length (n w ) of loops winding the honeycomb lattice as a function of lattice size 
L. A linear fit to the data obtained on lattices of size L = 32 up to L = 352 using 
(n w ) oc L D yields D = 1.37504(32) with xV D0F = 1.13, in excellent agree- 
ment with the predicted result. Note the extra digit of precision achieved here in 
comparison to the above estimates. 

In Fig. \\3\ we show the probability II l that one or more loops wind the lat- 
tice of size L at the critical point. The data show no finite-size effects. A one- 
parameter fit to the data with L ranging from L = 32 up to L = 352 gives 
I1 L = 0.51257(27) with xV D0F = i- 08 - Th is probability is slightly larger than 
the probability Il{Pc) = 0.50024(21) of finding a configuration with odd wind- 
ing number because 11^ also includes configurations with (nonzero) even winding 
number. 

To detail this further, we give in Table[5]the probability P w of finding a config- 
uration with winding number w measured on lattices of sizes L = 32 and L = 160 
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Figure 13: Probability IIl that one or more loops wind the honeycomb lattice of linear size L at 
the critical point. The most right symbol in the figure denotes the weighted average with error 
bars. 



at the critical point K 
with (N — 1) 



K c and also in the low-temperature phase at K = Kur 



K LT = [2 - (2 - AO 1 / 2 ] 



-1/2 



(53) 



defining the low-temperature branch of the O(iV) model on a honeycomb lattice 
lllql . At this temperature, the bond fugacity becomes unity for the Ising model so 
that each configuration carries the same weight according to the partition function 
(1421) . The loop model thus reduces to a purely geometric random model where 
a bond update is always accepted. By Eq. ( 1471 ), this temperature corresponds to 
vanishing inverse dual temperature $ — of the Ising model on the dual lattice 
where spins are oriented up or down at random. In this limit, the dual model 
becomes equivalent to random site percolation at the percolation threshold p c — |, 
and the polygons on the honeycomb lattice denote the boundaries of the occupied 
sites on the hexagonal lattice |33j]. For the two lattice sizes considered, we have 
not recorded any configuration with winding number w > 5. For the critical 
theory, we did not even observe a single configuration with w = 4. The data show 
no finite- size effects. The sum of the measured probabilities P w with w odd equals 
0.50032(89) for L = 32 and 0.49963(73) for L = 160 in the case of the critical 
theory, and 0.75016(41) for L = 32 and 0.74987(60) for L = 160 at /3 = 0. Since 
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Table 5: Probability of finding a closed loop configuration with winding number w on honeycomb 
lattices of size L = 32 and L = 160 at the critical point K c and at Kjjt where j3 = 0. 







K LT 


w 


= 32) 


P W (L = 160) 


P W (L = 32) 


P W (L = 160) 





0.48747(89) 


0.48802(74) 


0.15764(29) 


0.15836(29) 


1 


0.50029(88) 


0.49960(73) 


0.74579(34) 


0.74545(53) 


2 


0.01221(12) 


0.01236(14) 


0.09214(16) 


0.09171(36) 


3 


0.0000307(70) 


0.0000299(64) 


0.004369(70) 


0.004413(76) 


4 








0.0000588(64) 


0.000060(11) 



these numerical results are perfectly consistent with the fractions | and |, it is 
tempting to speculate that these are in fact exact results. 

As final observable, we consider the distribution function P„(r) of the end-to- 
end distance r introduced in Eq. ©. The top panel in Fig. [14] shows this distri- 
bution as a function of r for chains of length n = 615, 1230, and 1845 measured 
on the largest lattice considered, viz. L = 352. Each of the distributions are nor- 
malized according to Eq. (flOl) with d — 2. If our finite-size scaling conjecture © 
holds for the 0(1) loop gas, the data for various n should collapse onto a universal 
curve when n d / D P n with d = 2 and D = y is plotted as a function of the scaling 
variable r/n 1 ^ . This is indeed what we observe, see bottom panel in Fig. [14] 
The analyses of data measured on smaller lattices, typically with three equidistant 
chain lengths, give similar results. For SAWs it has been suggested [37] that the 
whole scaling function can be approximated by a single function 

V{t) = afexp (-bt 5 ) , (54) 

with parameters a and b. This scaling function slightly generalizes the form orig- 



inally proposed by Fisher II 1511 . where b = 1 was assumed. The exponent 5 is 



related to the fractal dimension D of the SAWs by the Fisher law IU5I1 

5 = -— . (55) 

We have cast this law in a form that allows generalization to arbitrary critical 
O(iV) loop gases, where it is recalled that for SAWs, which are described by the 
0(N 0) model, v = 1/aD with a = 1, but for arbitrary -2 < TV < 2, a ^ 1. 



u 

3 

data yields a surprisingly good approximation of the entire scaling function, see 



With 5 given the predicted value 5 = y for N = 1, a three-parameter fit to the 
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Table 6: Results of three-parameter fits to the data in Fig. [14] using the predicted form d54l l with 
6 = The last line in the table gives the weighted average of the estimates for L > 200, which 
are all consistent with a constant, i.e., L-independent value. 
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a 


b 





?9 'exact 


X 2 /dof 


160 


1.017(08) 


3.819(18) 


0.3970(73) 


1.058(19) 


1.80 


192 


0.992(11) 


3.798(24) 


0.377(10) 
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2.74 


224 


1.0070(94) 


3.820(19) 
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0.3793(56) 
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1.48 


288 


1.0022(50) 
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1.013(12) 


1.35 


320 


0.9950(47) 


3.790(10) 


0.3733(44) 


0.995(11) 


1.35 


352 


0.9973(51) 


3.814(10) 


0.3721(46) 


0.992(12) 


1.53 


oo 


0.9994(26) 


3.8096(55) 


0.3769(24) 


1.0050(64) 





Table [6] The resulting estimate for $ is in excellent agreement with the predicted 
value i? = §. Note that the normalization ( flOl) translates into the normalization 

poo 

l = tt d dtt d ~ l V n {t) (56) 
Jo 

of the scaling function. With the explicit form (1541) , this gives a relation between 
the three parameters a, b, and Our estimates for these parameters satisfy this 
relation within statistical errors, so that the fits effectively involve only two free 
parameters. 

6. Outlook and discussion 

The worm algorithm may potentially turn the loop gas approach to fluctuating 
fields on a lattice into a viable alternative for numerically simulating lattice field 
theories. Being an alternative, the algorithm calls for new estimators of physical 
observables, a few of which we have described here. Up to now, the loop gas 
approach has been widely adopted only in the de Gennes N — > limit of the 
O(N) model, which describes self-avoiding random walks. As demonstrated in 
this paper, concepts developed to describe such random walks can be generalized 
to arbitrary O(N) models. A next step in advancing the worm algorithm, which 
we hope to explore in a future publication, is to include gauge fields. Finally, 
although field theories that contain fermions can also be expanded in a strong- 
coupling series, and are in principle amenable to the worm algorithm, it is an 
open question whether the algorithm can deal with the sign problem. A first step 
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r / n yn 

Figure 14: Top panel: Distribution P n (r) (multiplied by a factor of 10 5 for convenience) as a 
function of the end-to-end distance r for chains of length n = 615, 1230, 1845 measured on the 
largest lattice considered, viz. L = 352. Bottom panel: Rescaled data shown in top panel. The 
curve through the data points is based on a three-parameter fit to the data using the predicted form 

da. 
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towards this problem, to which we also hope to return in a future publication, is 
to simulate the critical O(iV) model for N = — 1 on a honeycomb lattice, so that 
according to the HT representation (l42l) each polygon caries a minus sign. That is, 
in this particular representation, the 0(N = — 1) model exhibits the sign problem 
in its pristine form. 
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